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ABSTRACT 

We investigate the operation of the chemo-thermal instabihty in primordial star- 
forming clouds with a suite of three-dimensional, moving-mesh simulations. In line 
with previous studies, we find that the gas at the center of high-redshift minihaloes 
becomes chemo-thermally unstable as three-body reactions convert the atomic hydro- 
gen into a fully molecular gas. The competition between the increasing rate at which 
the gas cools and the increasing optical depth to H2 line emission creates a charac- 
teristic dip in the cooling time over the free-fall time on a scale of about 1000 au. As 
the effective equation of state drops from about 1.1 to well below unity, the free-fall 
time decreases to below the sound-crossing time, and density perturbations induced 
by transonic turbulence can grow. This allows the cloud to fragment on a scale of a 
few tens of au during the initial free-fall phase. In two of the nine haloes investigated, 
Jeans-unstable clumps condense out of the parent cloud, which will likely collapse 
before they are accreted by the primary clump. In the other haloes, fragmentation at 
such an early stage is less likely. However, given that previous simulations have shown 
that the infall velocity decreases substantially once the gas becomes rotationally sup- 
ported, the amount of time available for perturbations to develop may be much greater 
than is evident from the limited period of time simulated here. 

Key words: early Universe - cosmology: theory - galaxies: formation - galaxies: 
high-redshift - hydrodynamics - stars: formation 



1 INTRODUCTION 

The first stars in the Universe are believed to have formed 
only a few hundred million years after the Big Bang 



(Barkana h Loeb 2001 Bromm k, Larson 2004 Ciardi h 


Ferrara||2005| 


Glover||2005 Loeb||2010 Glover||2013 


). They 



heated and ionised the pristine intergalactic medium (IGM; 



Bromm et al. 2001; Schaerer 2002; Alvarez et al. 2006; Abel 
et al. 2007, Johnson et al ,2007, .Yoshida et al.,2007,. ; Whalen. 



et al.||2008 [Greif et al. p009| , and their supernova explo- 
sions enriched the primordial gas with the first heavy ele- 
ments (He ger et al.|[2003; Umeda h Nomoto||2003| [Wise 



Abel||2008| [Greif et al.||2010| ). They fundamentally altered 
the thermal and chemical state of the gas out of which the 
first galaxies formed, which in turn initiated the first self- 
sustaining cycle of star formation, feedback and chemical 
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enric hment ([Thoul h Weinb erg' 1996; Mac Low h Ferrara 



1999 


Madau et al.||2001| |0h & Haiman..2002| 


Ricotti et al. 


2002 


Wada & Venkatesan||2003| 


Dijkstra et a 


l.||2004| IRead 


et al. 


2006| |Wise & Abel|[2008| 


Greif et al. |2010|). Under- 



standing the formation and properties of the first stars is 
thus an important step towards a comprehensive picture of 
structure formation in the early Universe. 

The first theoretical studies of Population III star for- 
mation date back to the late 1960 's, in which molecular hy- 
drogen was recognised as a coolant of low-mass gas clumps 
that condensed out of the expanding Universe at high red- 
shift fSaslaw h Zipoy|19 67 ; Peebl es Dicke|1968 Hirasawa 
|T969 ; Matsuda et al.||1969 ; Takeda et al.||1969t . Subsequent 
studies used simplified one-zone models to account for the 
dynamics of collapsing gas clouds in the presence of radiative 
cooling ( [Yoneyam a 1972; Hutchins||1976| |Silk 



1977 



1983 



Carlberg 1981; Kashlinsky & Rees 1983; Palla et al. 1983" 
Carr et al.|1984,|Couchman ReeS|1986,|Uehara et al.|1996| 
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Tegmark et al. 1997), while the first one-dimensional calcula- 
tions of the simultaneous collapse of dark matter (DM) and 
gas were carried out in the context of the cold dark mat- 



ter (CDM) paradigm (Haiman et al.||1996| Omukai & Nishi 



1998|[Nakamura Umemura|1999 ). Three-dimensional sim- 
ulations of Population III star formation had to await im- 
provements in numerical simulation techniques, and were 



not possible until the late 1990 's ( [Abel et al, 
2002||Bromm et al.||1999| |2002| ) 



1998 2000 



Based on these pioneering studies, a now widely ac- 
cepted 'standard model' of Population III star formation has 
emerged. In this picture, the first bound structures capable 
of hosting star formation are considered to be DM 'mini- 
haloes' with virial masses Mvir ^ 10^ M©, which collapse at 
redshifts z > 20. At the center of these haloes, molecular 
hydrogen forms due to the non-negligible free electron frac- 
tion left over after recombination. The ro-vibrational states 
of H2 are excited by collisions with atomic hydrogen and 
helium, which decay radiatively and cool the gas. The min- 
imum temperature of about 200 K is reached at a density 
of ^ 10^ cm~^, where the level populations of H2 transition 
into local thermal equilib rium (LTE; |Abel et al.|1998[|2002| 
Bromm et al.|1999 2002). Shortly thereafter, the central gas 
cloud becomes Jeans-unstable and collapses in its own right. 

At densities nu > lO^cm"^, where nu is the num- 
ber density of hydrogen nuclei, three-body reactions convert 



the remaining atomic hydrogen into H2 ( 


Abel et al. 


2002 


Bromm & Loeb|2004| Yoshida et al.|2006t 


Turk et al.|2009). 



During this phase, the gas may become chemo-thermally 
unstable, since the rapidly increasing H2 fraction results 
in increased cooling, which further accelerates the collapse 
and the rate at which H2 is formed ( Sabano Yoshii|1977 



Silk 1983 Omukai & Yoshii 2003). This runaway process 



is counteracted by the release of the binding energy of the 
formed H2 , and the increasing optical depth of the gas to H2 
line emission. As the molecular hydrogen fraction saturates, 
collision-induced emission becomes important, followed by 
collisional dissociation of H2 (Ripamonti & Abelf2004) . The 
final, nearly adiabtic contraction of the cloud ends with the 
formation of a protostar with a mass of c^i 0.01 Mq ( [Omukai] 
k Nishi||1998| [Yoshida et al.|[2b08( ). 

An inherent problem of self- gravitating collapse simula- 
tions is the continually decreasing timescale on which the gas 
evolves, which limits the integration time to of order a few 
dynamical times of the Jeans-unstable cloud. This problem 
may be circumvented by employing so-called sink particles, 
but comes at the cost of introducing artificial boundary con- 



ditions ( Bate et al.|1995|[Krumholz et al.|2004[[Jappsen et al.[ 
2005 Federrath et al.[[2010 ). Recent studies that employed 
sink particles found that the gas becomes rotationally sup- 
ported in a Keplerian disc. The disc becomes gravitationally 
unstable due to the high mass accretion rate from the en- 
velope onto the disc and the efficient cooling of the disc by 
H2 line emission, and fragments into a handful of protostars 
(Clark et al. 2008, 2011a, b; Stacy et al. 2010). A follow-up 
study that investigated a sample of minihaloes found that 
the star- forming clouds in them fragmented in a similar fash- 
ion ( jGreif et al.(20TT| . 

Simulations that avoid sink particles and instead di- 
rectly resolve the collapse of the gas to protostellar densities 
have become possible, but their high computational limits 
the integration time to ~ 10 yr ( Greif et al.|2012 ). In analogy 



to previous simulations with sink particles, this study finds 
that the circumstellar disc fragments into a small system of 
protostars. Strong torques exerted by the protostars and the 
disc result in the migration of about half of the protostars 
to the center of the cloud in a free-fall time, where they 
merge with the primary protostar. A smaller fraction of the 
protostars migrates to higher orbits and survives until the 
end of the simulation. Two-dimensional simulations with a 
piecewise polytropic equation of state show similar patterns 
of migration and merging fVorobyov et al. [[2013 ). 

The susceptibility of the gas to fragmentation is a key 
ingredient towards understanding the initial mass function 
of the first stars. An early termination of fragmentation may 
result in the formation of a single, massive Population III 
star, while continued fragmentation may lead to the for- 
mation of a rich cluster of less massive stars. Most of the 
mass that is accreted before the gas is dispersed by ionis- 
ing radiation resides betw een 10 and 1000 au (^Tan McKee 



2004)[McKee Tan|2008| [Hosokawa et al.|2011| [Stacy et aL 



2012), implying that it is crucial to understand the thermal 
evolution of the gas on these scales. In the present study, 
we use a suite of three-dimensional, moving mesh simula- 
tions to investigate the collapse of the gas on these scales in 
detail. Compared to previous work, we employ higher reso- 
lution and investigate a larger sample of haloes. This allows 
us to assess the properties of primordial gas clouds with 
unprecedented accuracy. In addition, we perform the most 
comprehensive resolution study of primordial star formation 
to date, using up to three billion resolution elements. 

The structure of our work is as follows. In Section 2, we 
describe the numerical methodology, the setup of the simu- 
lations, and the chemistry and cooling network. In Section 3, 
we present the results of a one-zone calculation that high- 
lights the various physical processes that operate on the rele- 
vant scales, discuss the thermal and gravitational stability of 
the gas in a representative halo in the three-dimensional sim- 
ulations, followed by an analysis of the star-forming clouds 
in all haloes. We then discuss the results of the resolution 
study, and in Section 4 summarise our main findings and 
draw conclusions. All distances are quoted in proper units, 
unless noted otherwise. 



2 SIMULATIONS 

The three-dimensional, cosmological hydrodynamics equa- 
tions for a mixed fluid consisting of coUisionless DM and 
collisional gas are solved using the methodology employed 
in the moving-mesh code AREPO (Springel 2010). In the fol- 
lowing, we describe in detail the initialisation and setup of 
the simulations. 



2.1 Dark matter simulations 

The background Universe used to initialise the simulations 
is a standard A cold dark matter (ACDM) cosmology with 
parameters based on the PI/MAP observations (e.g. K omatsu 
[et al.|2009t . We use a matter density of Qm = 1 



-1}a = 0.27, 

baryon density = 0.046, Hubble parameter h — 0.7, spec- 
tral index rig = 0.96, and normalisation ag = 0.81. The mat- 
ter power spectrum is evolved forward in time until 2; = 99, 
and we use the Zel'dovich approximation to determine the 
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initial displacements of the particles that represent the DM 
distribution function from a cubical lattice. We initialise 
nine boxes with a side length of 1 Mpc (comoving) and dif- 
fering realisations of the power spectrum to investigate how 
our results are affected by cosmic variance. We employ 512^ 
particles of mass ^ 272 , and use a gravitational soft- 
ening length of ^ 98 pc (comoving), which corresponds to 
5% of the initial mean inter-particle separation. The simu- 
lations are evolved until the first halo grows to a virial mass 
of Mvir = 5 X 10^ M0, which is evaluated by an on-the-fly 
friends-of- friends algorithm ( Springel et al.|2QQlJ . 



Table 1. Halo properties 



2.2 Resimulations 

Once the halo designated for further refinement has been 
located, the simulations are centred on this halo and reini- 
tialised with higher resolution. For this purpose the particles 
in the target halo and a sufficiently large boundary region 
around it are traced back to their initial positions, which 
yields the Lagrangian volume out of which the halo formed. 
In this region, the resolution is increased by a factor of 16^ 
with respect to the original simulation, and augmented with 
additional small-scale power. Each DM particle is replaced 
by a less massive DM particle and gas cell with a displace- 
ment corresponding to half of the initial mean inter-particle 
separation, and a mass ratio of Mgas/Mdm = ^h/{^m — ^h)- 
Outside of the target region, the resolution is gradually de- 
creased to reduce the total number of resolution elements, 
while the accuracy of the gravitational tidal field that in- 
fluences the formation of the halo is preserved. The DM 
particle and gas cell masses in the high-resolution region are 
~ 0.05 and c^i 0.01 M©, respectively, corresponding to nearly 



100 times higher resolution compared to Greif et al. (2011 



2012). The gravitational softening length is set to ^ 6pc 
(comoving) . 

The cosmological resimulations are run until the min- 
imum cell size is of order 10~^ of the side length of the 
box. In the present version of AREPO, the number of exact 
integer operations used to construct a valid Voronoi tesse- 
lation increases substantially when the dynamic range ex- 
ceeds the above value, which decreases the performance of 
the code. Nearly simultaneously, the maximum resolution 
of the Peano-Hilbert curve is reached and the cells can no 
longer be distributed efficiently among tasks. Both of these 
limitations will be addressed in a future version of the code. 
Here, we resort to extracting the central 1 kpc (comoving) of 
the box and reinitialising the simulations at the final output 
time. Once the density exceeds nu = 10^ cm~^, we perform 
a second extraction and cut out the central 1 pc, discard 
the DM particles, and evolve the simulations to a density 
of riH = 10^^ cm~^. The intermediate resimulations are nec- 
essary to achieve a maximum dynamic range in the final 
resimulations, which simplifies the analysis of the data. In 
both cases, the sound crossing time through the box is much 
longer than the free-fall time of the dense gas at the center, 
such that perturbations from the edge of the box do not 
affect the central cloud. 

In terms of computational efficiency, AREPO has been 
further improved with respect to the version employed in 
Greif et al. (2012). One of the most important new features 
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Redshift of collapse, virial mass, virial radius and spin param- 
eter of the first minihaloes that form in the simulations. 

ancing by taking into account the wall clock time spent on 
the various levels in the time step hierarchy. Furthermore, 
the gravitational oct-tree is now reconstructed at each time 
step, which prevents a dramatic increase in the wall clock 
time spent on the tree walk when the sizes of the nodes 
are adjusted dynamically to encompass all particles that are 
assigned to them. Finally, a separate neighbour tree is con- 
structed to speed up the neighbour search used in particular 
for the construction of the Voronoi tesselation. 



2.3 Refinement 

The so-called Truelove criterion states that at least four cells 
per Jeans length are necessary to resolve gravitational insta- 
bility ( Truelove et al.| p^998). However, this is generally not 
sufficient to resolve the turbulent cascade that is driven by 



gravitational collapse. For example, Federrath et al. (2011) 
showed that at least 32 cells per Jeans length are necessary 
to model the amplification of magnetic fields by the turbu- 
lent dynamo. In Turk et al.| (.2012), it was shown that the 
employed resolution affects the chemical, thermal and dy- 
namical evolution of primordial gas clouds. In particular, 
the amount of rotational support increases as the resolution 
is decreased, which might affect the ability of the gas to 
fragment. Although the resolution requirements in a fixed 
grid simulation may not be the same as in a moving-mesh 
simulation, these results show that it is desirable to use as 
many cells per Jeans length as possible. 

[Turk et al.| ( |2010| ) further found that it is advantageous 
to evaluate the Jeans length using the minimum tempera- 
ture of the gas, which ensures that shock-heated regions are 
equally well resolved as cold regions. To facilitate a com- 
parison with previous studies, we therefore also employ the 
refinement criterion proposed by |Turk et al.| ( |2010| ), and eval- 
uate the Jeans length at Tmin = 200 K, which yields: 

/ nx. \ -1/2 
Aj :^0.6pc(^: 



. 10^ cm- 



(1) 



is a domain decomposition optimised for simulations with a 
high dynamic range, which allows close to optimal work bal- 



Using this definition, we refine cells if h > Xj/Nj, where 
h = {SV/Any^^ is the approximate radius of a cell, V its 
volume, and A^j the desired number of cells per Jeans length. 
For A^j = 32 and nu > 10^ cm~^, the Jeans mass is resolved 
by 10 times more cells compared to [ Greif et al.| p012). In 
Section 3.3, we present a resolution study in which we sys- 
tematically increase A^j from 8 to 128. Due to the excessive 
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Figure 1. From top left to bottom right: temperature, H2 frac- 
tion, effective equation of state, and cooling time over free-fall 
time versus number density of hydrogen nuclei in a one-zone cal- 
culation. The various profiles show the evolution for a constant H2 
fraction (blue), indefinitely increasing H2 fraction (green), limited 
H2 fraction (red), three-body formation heating (cyan), and an 
increasing optical depth to line emission (magenta). The dashed 
grey lines denote an effective equation of state of 1.1 and a cool- 
ing time equal to the free-fall time. If the H2 fraction is allowed 
to increase indefinitely, the chemo-thermal instability results in 
an asymptotic 7eff = 1/2. However, the decreasing H2 fraction, 
three-body formation heating, and the optical depth of the gas 
reduce the effectiveness of the instability, such that the cooling 
time over the free-fall time drops only slightly below the value 
expected for a constant H2 fraction. 



computational cost of simulations with A^j = 128, we use 64 
cells per Jeans length for the main simulation runs discussed 
in Sections 3.1 and 3.2. 

In addition to the refinement based on the Jeans length, 
we refine cells if their mass increases to more than twice their 
initial mass. While AREPO maintains near constant mass res- 
olution by advecting the cells in a Lagrangian fashion, the 
Eulerian mass flux between cells may lead to a gradual drift 
away from the initial masses of the cells. 



2.4 Chemistry and cooling 

The chemistry and cooling network employed in the present 
study is nearly identical to the one described in | Greif et al.| 
([2012^. The species evolved next to the internal energy of the 
gas are H, H+, H", H+, H2, He, He+, He++, D, D+, HP, and 
free electrons (Glover & Jappsen 2007 Clark et al. 2011). 



We use the publicly available SUNDIALS CVODE package 
to solve the coupled chemical and thermal rate equations 
(Hindmarsch et al. 2005). The abundances are computed 
with a relative accuracy of 0.1%, which is necessary for an 
accurate integration at high densities, where three-body re- 
actions become important. The minimum abundance that is 
computed reliably is yx = 10~^°, where yx denotes the num- 
ber density of species X with respect to the number density 
of hydrogen nuclei. This threshold prevents the solver from 
computing accurate abundances when it is not necessary. 
The most important low-temperature coolant in pri- 



Figure 2. From top left to bottom right: number density of hy- 
drogen nuclei, H2 fraction, photon escape fraction, and cooling 
time over free-fall time versus radius. The dashed grey line de- 
notes a cooling time equal to the free-fall time. The H2 fraction 
rapidly increases below ~ 1000 au due to three-body reactions, 
which results in the cooling time falling slightly below the free- 
fall time. As the optical depth of the gas to H2 line emission 
increases, the cooling rate decreases again, and the cooling time 
remains nearly equal to the free-fall time down to scales of a few 
tens of au. 



mordial gas is molecular hydrogen, which forms via associa- 
tive detachment of H and H~ , following radiative association 



of H and e ( 


Saslaw k Zipoy||l967 


Abel et al.||l997 


Gain 


& Palla|1998) 


. The ro- vibrational states of H2 are collision- 



ally excited and decay radiatively, allowing the gas to cool 
to a minimum temperature of T :^ 200 K. At low densities, 
we use the excitation rates for collisions with neutral hydro- 
gen atoms given by Wrathmall & Flower ( 2007| , assuming 



an ortho-para ratio of 3:1. At high densities, when the en- 
ergy states are populated according to LTE, we use the rates 



given by Le Bourlot et al. (1999). We neglect collisions with 
electrons, protons, and neutral helium atoms, which may 
provide a no n- negligible amount of cooling at the relevant 
densities and temperatures ( Glover &: Abel||2008t . 

Another low-temperature coolant in primordial gas is 
hydrogen deuteride (HD), which may become important at 
temperatures T < 200 K due to its permanent dipole mo- 
ment (Stancil et al.||1998| [Flower ~al. 2000 j Nakamura 



Umemura||2002| [Johnson "fc Bromm||2006r|R]pamonti||2007[ 



McGreer Bryan[[2008[ [Greif et al.[[201l| . In minihalos, 
HD forms mainly via associative detachment of H2 and D^, 
where the abundance of D^ is set mainly by recombina- 
tions and charge transfer with H. The cooling rate of the 
gas by HD ro-vibrational transitions is given by jLipovka[ 
et al.[ ( [2005t . 

At densities nu > 10^ cm""^, molecular hydrogen forms 
via three-body reactions, which heats the gas due to the 



release of the molecular binding energ y (Palla et al 



Bromm et al.[[1999 



[Abel et al.|[2002l [Bromm et ~ 



1983 



2002 



Yoshida et al.[[2006 ). The corresponding reaction rates are 
highly uncertain, which is reflected by the substantial varia- 
tion of the properties of primordial gas clouds when the dif- 
ferent rates are used ( ^Turk et al.^^2011] . We here adopt the 
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read datasets in C and distribute the particle load to the 
processing cores. The operations required for the various 
tasks, e.g. the image generation or the computation of ra- 
dial profiles, are then performed on the particles residing on 
each core, allowing near-optimal work balance. Upon com- 
pletion, an intermediate dataset is written that is pipelined 
to a Python program for further analysis. The combined 
workflow of SATOR is controlled by a script that reads in a 
minimal set of parameters from the command line, which 
specifies the mode of operation and the target dataset. The 
script may be run interactively, or submitted to a queuing 
system for very large datasets and time-consuming calcula- 
tions. On a Sandy bridge computing cluster, the runtime of 
SATOR to analyze our largest dataset with ~ 3 x 10^ cells 
using 512 MPI tasks is of order 10 s, which is a substantial 



improvement compared to Greif et al. (2011 2012) 



Figure 3. From top left to bottom right: temperature, effec- 
tive equation of state, free-fall time over sound-crossing time, and 
RMS density contrast versus radius. The dashed grey lines denote 
an effective equation of state of 1.1, a free-fall time equal to the 
sound-crossing time, and an RMS density contrast of unity. As the 
cooling time decreases, the effective equation of state of the gas 
drops from 7eff > 1.1 to ~ 0.8, resulting in a slight drop in tem- 
perature from ~ 1000 K to ~ 800 K below ~ 1000 au. As a result, 
the free-fall time falls below the sound-crossing time and pertur- 
bations seeded by the transonic turbulence within the cloud can 
grow. This is evident from the bottom right panel, which shows 
that the RMS density contrast increases from about as — 0.6 on 
a scale of ~ 10^ au to well above unity on a scale of ~ 100 au. 



intermediate rate for three-body H2 formation among those 
(2011|, which is taken fromlPalla et al.l (11983). 



of [Turk et ah 

For flu ^ 10^^ cm"*^, column densities are high enough that 
the gas becomes optically thick to H2 hne emission, and we 
use the Sobolev approximation to determine the escape frac- 
tion of photons emitted by H2 lines ( jYoshida et al.|2Q06J . 

In contrast to Greif et al. ( 2012 ), we do not switch to an 
equilibrium chemistry solver at a density of nu = 10^^ cm~^, 
since we only follow the collapse up to a density of nu — 
10^^ cm"*^. However, we include a number of performance 
optimisations that decrease the computational cost of the 
non-equilibrium solver by a factor of a few. In particular, we 
improve cache utilisation by minimising the time spent on 
the retrieval of reaction rates, which are accessed repeatedly 
as the abundances are evolved. Finally, our main simulations 
do not include species transfer between cells, which has only 
recently been implemented in AREPO. However, we have per- 
formed a test simulation in which we follow the high-density 
evolution of the gas in MHl using species transfer, and found 
that it only marginally affects our results. 



2.5 Analysis 

The number of resolution elements employed in the main 
simulations exceeds 10^, which is an increase by more than 
two orders of magnitude compared to the simulations pre- 



sented in Greif et al. (2011 2012). Such large datasets also 
require more sophisticated tools for an efficient analysis. 
To this end, we have written the parallel software package 
SATOR, which uses the message passing interface (MPI) to 



3 RESULTS 

3.1 Chemo-thermal instability 

3.1.1 One-zone model 

To qualitatively understand the influence of the various 
physical processes on the relevant scales, we perform a one- 
zone collapse calculation that employs a simplified chemistry 
and cooling model. In a roughly isothermally contracting 
cloud, the rate at which the core density increases is given 



by: 

nn = nn/tff, 



(2) 



where is = a/Stt/ (32Gp), G is Newton's constant, p — 
munu/X is the mass density, mu is the mass of the hy- 
drogen atom, and X = 3/4 is the cosmological mass frac- 
tion of hydrogen. The volumetric heating rate of the gas in 



erg s cm 



is given by: 
Ayu2f esc i^- 



(^) +^^"^^2 



(3) 



where 7 is the adiabatic index of the gas, 2/H2 the molec- 
ular hydrogen fraction, A = 10~^^ erg cm~^ s~^ and a — A 
approximate the H2 Hne cooling function used in the full 
three-dimensional simulations, B — 4.48 eV is the binding 
energy of the H2 molecule, and /esc the photon escape frac- 
tion for H2 Hne emission. The first term on the right-hand 
side of the equation describes adiabatic heating, the sec- 
ond term H2 Hne emission, and the third term three-body 
formation heating. We have verified that all other thermal 
processes are comparatively sub-dominant at the relevant 
densities and temperatures. We use an adiabatic index of 
7 = 5/3 for a monatomic gas, even though the gas consists 
of a mix of atomic and molecular hydrogen. We also neglect 
the change in the mean molecular weight /x, which counter- 
acts the effect of the varying 7 in the relation between the 
energy density and temperature. The net effect of neglecting 
both processes is thus small. The photon escape fraction is 
given by: 



1.45 



3,1.45 _^ 0.45 



(4) 



for X ^ 1 and /esc = 1 for x < 1, where x — tih/'/T'H, thresh 
and nn, thresh = 4 x 10^ cm~^. This formula is based on the 
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Figure 4. From top left to bottom right: number density of hydrogen nuclei, temperature, H2 fraction, and free-fall time over sound- 
crossing time in a cube of side length 200 au, weighted with the mass and the square of the density of the cells traversed along the line of 
sight. As the chemo-thermal instability develops and the cooling time falls below the free-fall time, a secondary clump condenses out of 
the parent cloud that is visible on the left-hand side of the panels. The gravitational instability of the clump is evident from the bottom 
right panel, which shows that the free-fall time has dropped to well below the sound-crossing time. 



fit provided by Ripamonti & Abel (2004), but has the ad- 



vantage of a continuous derivative at x = 1. Finally, the rate 
at which H2 molecules are formed is given by 



(5) 



where C = 5.5 x 10 



We have again verified that all 



other rate equations that affect the chemical composition of 
the gas are comparatively small at the relevant densities and 
temperatures. We initialise the calculations at a density of 
riH = 10^ cm~^, using a temperature of T = 580 K and an H2 
fraction of = 10~^. These initial conditions roughly cor- 
respond to the expected values at the initial density, and in 



addition yield a smooth transition into the regime at which 
the chemo-thermal instability operates. 

In Figure [1] we show the results of the one-zone calcula- 
tion as we systematically increase the level of physical detail. 
The blue line shows the evolution for a constant H2 fraction 
and a constant escape fraction of unity. In this case, the ef- 
fective equation of state remains close to 7eff = 7/6, since 

Q 1/2 rn 1 / m 

tcooi — tff, which yields T~ oc and thus T (x p ^ . The 

green line shows the collapse assuming that the H2 fraction 
may increase indefinitely, the three-body formation rate is 
independent of the H2 fraction, and that three-body forma- 
tion heating as well as the influence of the decreasing es- 
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Figure 5. From top left to bottom right: temperature, escape 
fraction, cooling time over free-fall time, and free-fall time over 
sound-crossing time versus number density of hydrogen nuclei. 
The logarithm of the mass per bin over the total mass in the 
box is color-coded from blue (lowest) to red (highest), and the 
solid lines show the mass-weighted average value versus density. 
The grey dashed lines denote a cooling time equal to the free- 
fall time, and a free-fall time equal to the sound-crossing time. 
The temperature distribution widens significantly above nn ~ 
10^° cm~^ as the cooling time falls below the free-fall time. The 
photon escape fraction at a given density fiuctuates by up to two 
orders of magnitude, such that some regions of the cloud cool 
very efficiently and may become Jeans-unstable. This is reffected 
by the substantial amounts of gas with a free-fall time below the 
sound-crossing time. 



cape fraction may be neglected. In this case, an asymptotic 
7eff = 1/2 is reached, which is evident from the following 
calculation. 

The cooling time due to H2 line emission in the optically 
thin limit is given by: 

tcool oc ^HsT"^, (6) 

and the derivative of the H2 fraction with respect to nn is 
given by: 

dyu2 _ d2/H2 dt 



dnu dt dnu 

where we have assumed T oc p^^^^"^ 
and insertion into equation (6) yields: 



(7) 

Integration over riHa 
(8) 



, 3/2-27eff 
t'cool oc TT'ii^ 1 

and the cooling time over the free-fall time is thus given by: 



^cool 



l-27eff 



(9) 



A constant ratio of the cooling time to the free-fall time is 
therefore only possible for 7eff = 1/2, which agrees with the 
results shown in Figure [l] 

Under more realistic circumstances, this asymptotic 
value is never achieved, since only a finite amount of H2 
may be formed. In addition, the rate at which H2 forms de- 
creases as the H2 abundance increases. The red line in Fig- 
ure [1] denotes the case in which these considerations have 
been taken into account. As expected, the effective equation 



of state departs from the asymptotic 7eff = 0.5 as the H2 
abundance increases to of order unity, and quickly returns 
to 7eff = 7/6. When three-body formation heating and the 
increasing optical depth to H2 line emission are taken into 
account (cyan and magenta lines, respectively), the effec- 
tiveness of the chemo-thermal instability is further reduced, 
and the temperature, effective equation of state, and ratio 
of the cooling time to the free-fall time are significantly less 
affected. Nevertheless, an unstable regime develops at a den- 
sity of > riH — 10^ cm""^, where 7eff drops to 0.9. 



3.1.2 Simulations 

The properties of the gas in the three-dimensional simu- 
lations are similar, but not entirely equal to those of the 
one-zone calculations. In Figure |2] we show the number den- 
sity of hydrogen nuclei, H2 fraction, photon escape fraction, 
and cooling time over free-fall time versus radius r at the 
last output time of the high-density simulations in MHl. 
The profiles are determined by using a mass- weighted aver- 
age of the cells that contribute to the radial bins. The top 
right panel shows that the H2 fraction increases rapidly on 
a scale of 1000 au as three-body reactions become impor- 
tant. This results in the cooling time decreasing to about 
half of the free-fall time, which is evident from the dip at 
r 1000 au. On a scale of 500 au, the optical depth of 
the gas to H2 line emission increases, and the cooling time 
again approaches the free-fall time. 

In Figure[3] we show the temperature, effective equation 
of state, free-fall time over sound-crossing time, and root- 
mean-square (RMS) density contrast versus radius. The lat- 
ter is given by: 



E 



(10) 



where the sum extends over all cells contributing to a radial 
bin, i denotes the cell under consideration, rrii its mass, pi 
its density, Mtot the total mass per bin, and p the mass- 
weighted average density of the bin. Following the decrease 
of the cooling time below the free-fall time, the top left panel 
in Figure |3] shows that the temperature of the gas briefly 
decreases from 1000 K to ^ 800 K, which is reflected in 
the drop of the effective equation of state from 7eff > 1.1 
to :^ 0.8. For 7eff < 7crit — 1-1, it is known that density 
perturbations may be amplified as the gas further contracts 
( Hanawa Matsumoto|[2(300| ) . The growth rate of the per- 
turbations depends on the mode of the perturbation and 
the exact value of the effective equation of state. Generally, 
the largest modes grow most rapidly, but the achieved den- 
sity contrast also depends on the initial amplitude of the 
perturbations. The transonic turbulence that permeates the 
cloud seeds the strongest perturbations on small scales, but 
also washes them out most quickly on these scales. In ad- 
dition, the growth of small-scale perturbations below the 
Jeans length is suppressed by pressure gradients. The most 
rapidly growing perturbations therefore reside on scales just 
above the Jeans length. 

Once the cooling time falls below the free-fall time, the 
sound-crossing time decreases more rapidly than the free- 
fall time. Perturbations can thus develop more readily and 
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Figure 6. Number density of hydrogen nuclei in cubes of side length 200 au, centred on the densest cell in each halo. The number density 
is weighted with the mass and the square of the density of the cells traversed along the line of sight. Overall, the clouds are centrally 
concentrated, but display very irregular morphologies with filaments and knots indicative of turbulence. In MH5 and MH7, elongated 
clouds have formed, while in MH6 a prominent disc with little substructure is present. In MH5 and MH7, the collapse of a secondary 
clump has progressed far enough that it is clearly discernible from the first. Somewhat less pronounced is the sub-clump that has formed 
in MHl. 



become apparent within the free-fall time of the primary 
clump. A minimum in the free-fall time over the sound- 
crossing time develops on scales significantly below the mini- 
mum in the cooling time over the free-fall time, since a signif- 
icant change in the sound speed requires of order a cooling 
time, which drops to a minimum of about tff/2 (see Fig- 
ure |2]). This is reflected in the minimum of the free-fall time 
over the sound-crossing time being located at a few tens of 
au, while the minimum in the cooling time over the free-fall 
time is located at :^ 500 au. The elevated density contrast 
that develops as perturbations grow is evident from the bot- 
tom right panel, which shows that the RMS density contrast 



increases from about cr^ :^ 0.6 on a scale of ^ 10^ au to well 
above unity on a scale of :^ 100 au. 

The tentative formation of a secondary clump during 
this stage is evident from Figure [4j which shows the number 
density of hydrogen nuclei, temperature, H2 fraction, and 
free-fall time over sound-crossing time in the central 200 au 
of the box. The plane of view is perpendicular to the net an- 
gular momentum vector of the gas. A fully molecular clump 
with a temperature of ^ TOOK has formed that is visible 
on the left-hand side of the panels. The potential gravita- 
tional instability of the clump is apparent from the bottom 
right panel, which shows that the free-fall time has dropped 
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Figure 7. From top left to bottom right: number density of hydrogen nuclei, enclosed gas mass, temperature, H2 fraction, radial velocity 
over sound speed, rotation velocity over sound speed, rotation velocity over Keplerian velocity, turbulent Mach number, cooling time 
over free-fall time, effective equation of state, free-fall time over sound-crossing time, and RMS density contrast versus radius. The solid 
black line in the top left panel denotes an r"-^ -^ density profile. The grey dashed lines denote a cooling time equal to the free-fall time, 
an effective equation of state of 1.1, a free-fall time equal to the sound-crossing time, and an RMS density contrast of unity. A detailed 
discussion of this figure is provided in Section 3.2. 
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Figure 8. Free-fall time over sound-crossing time in cubes of side length 200 au, centred on the densest cell in each halo and weighted 
with the mass and the square of the density of the cells traversed along the line of sight. Parcels of gas that have a tendency to become 
Jeans-unstable are visible in MHl, MH5, MH7 and MH9. A comparison with Figure |6] reveals that they are are correlated with the 
isolated, dense clumps visible in this figure. 



to well below the sound-crossing time. In Section 3.3, we 
investigate whether this clump has acquired enough mass 
to become Jeans-unstable in its own right, and whether it 
collapses before it is accreted by the primary clump. 

The growth of perturbations enabled by the chemo- 
thermal instability is further evident from Figure [5] which 
shows the temperature, escape fraction, cooling time over 
free-fall time, and free-fall time over sound-crossing time 
versus number density of hydrogen nuclei. The logarithm 
of the mass per bin over the total mass in the box is color- 
coded from blue (lowest) to red (highest), and the solid lines 
show the mass- weighted average value versus density. Above 
a density of nu ^ 10^°cm~^, the temperature distribution 
widens and some parcels of gas cool to as low as ^ 200 K, 



while others are heated to over 2000 K. The photon escape 
fraction at a given density fluctuates by up to two orders 
of magnitude, such that some regions of the cloud cool very 
efficiently and may become Jeans-unstable. This is reflected 
by the substantial amounts of gas with a free-fall time be- 
low the sound-crossing time. Another noteworthy feature is 
the shock- heating of a parcel of gas to over 3000 K, which is 
apparent from the peak in temperature between nu 10^ 
and ^ 10^°cm~^. The elevated temperature results in col- 
lisional dissociation of H2 down to a level of 2/H2 — 10~^. 
Similar shock- heated parcels of gas were also found in |Turk| 



et al. (2010) 
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divided by the total mass Mtot 
Keplerian velocity is given by vi^ep — GM^^ 
Menc is the enclosed mass. The turbulent velocity is given 
by 



and radius of t he b in. The 
c/r, where 



1.0 1.5 

log r [au] 

Figure 9. Enclosed mass over Bonnor-Ebert mass (top panel), 
and free-fall time over accretion time (bottom panel) versus ra- 
dius centred on the secondary clumps that have formed in MHl, 
MH5, MH7 and MH9. The vertical lines denote the separation of 
the clumps from the center of the cloud. Only MH5 and MH7 are 
likely to fragment during the initial free-fall phase, since the ra- 
dius at which the enclosed mass exceeds the Bonnor-Ebert mass 
and the free-fall time drops to below the sound-crossing is signif- 
icantly smaller than the separation of the clumps from the center 
of the cloud. This is not the case in MHl and MH9, where the 
primary and secondary clumps are nearly indistinguishable. 



3.2 Halo sample 

In Figure [6] we show the number density of hydrogen nu- 
clei in a face-on views of the central 200 au in the various 
haloes investigated. Most clouds display a centrally concen- 
trated, irregular morphology with a multitude of filaments 
and knots, while MH4 and MH6 have a more regular mor- 
phology and less substructure. In MH5 and MH7, an elon- 
gated cloud has formed, while in MH6 a prominent disc with 
little substructure is present. In MH5 and MH7, the collapse 
of a secondary clump has progressed far enough that it is 
clearly discernible from the first. Somewhat less pronounced 
is the sub-clump that has formed in MHl. 

Figure [7| shows radially averaged profiles of various 
physical quantities. From top left to bottom right, the pan- 
els show the number density of hydrogen nuclei, enclosed gas 
mass, temperature, H2 fraction, radial velocity over sound 
speed, rotation velocity over sound speed, rotation veloc- 
ity over Keplerian velocity, turbulent Mach number, cooling 
time over free-fall time, effective equation of state, free-fall 
time over sound-crossing time, and RMS density contrast 
versus radius. The rotation velocity vector Vrot is given by 
the sum of the angular momenta of the contributing cells 



'^turb 



m^ {Vi - Vrad - Vrot) 



(11) 



where Vi is the velocity vector of a cell, and Vrad the radial 
velocity vector of the bin. The turbulent Mach number A4 
is given hy M = '^^turb/cs, where Cs is the sound speed. 

The top left panel in Figure [7| shows that the gas clouds 
have an approximate r~^'^ profile, which is expected for 
7eff ^ 1.1 ( |Larson|[T969 ). In MH5 and MH7, bumps in the 
density profile at a distance of a few tens of au from the cen- 
ter indicate that secondary clumps have begun to form. The 
enclosed mass is less sensitive to fluctuations in the density 
and does not show any prominent features. The temperature 
rises to about 1000 K on a scale of ^ 1000 au, followed by 
a slight drop to 800 K on a scale of 10 au, and an increase 
to over 1000 K on scales below :^ 1 au. In MH6, HD cooling 
becomes important and prolonges the initial cooling phase 
to riH — lO^cm"^, after which the temperature increases 
more sharply compared to haloes in which HD cooling does 
not become important. This leads to a slower collapse of 
the cloud during which turbulent motions decay and a pro- 
nounced disc forms (| Clark et al. 2011 ). The molecular hydro- 
gen fraction shows little scatter among the individual haloes 
and increases from 2/H2 — 10""^ at ^ 1000 au to fully molec- 
ular with yii2 — 0.5 at a few tens of au. The radial velocity 
peaks on a scale of a few tens of au as the gas cools slightly 
and loses some pressure support. 

All clouds are substantially rotationally supported. The 
average rotation velocity is of order the sound speed and 
comparable to the radial velocity. Most haloes display a peak 
in the respective ratio on a scale of a few tens of au, which 
is however not present in the profiles showing the rotation 
velocity over the Keplerian velocity. In MH6, the cooling 
time briefly increases to nearly an order of magnitude above 
the free-fall time, and the prominent disc visible in Figure |6] 
forms. This is reflected by an increase of the rotation velocity 
over the Keplerian velocity to 0.6 on a scale of ^ 1000 au. 
The respective value fluctuates by about a factor of two in 
the course of the collapse of a cloud, and the mean among 
the haloes is :^ 1/3. The turbulent Mach number is generally 
close to unity, implying that the clouds are in nearly equal 
parts supported by thermal pressure, rotation and turbu- 
lence. Below ^ 100 au, the temperature drops slightly and 
the turbulence transitions into the supersonic regime with 

~ 2. The turbulent support is also apparent from the 
filamentary appearance of the gas shown in Figure |6] 

Following the rapid increase of the H2 fraction on a scale 
of :^ 1000 au, the cooling rate increases and the cooling time 
drops slightly below the free-fall time. The effective equa- 
tion of state of the gas decreases from 7eff — 1.1 to ^ 0.8 on 
a scale of 100 au, and the free-fall time drops below the 
sound-crossing time, with MHl, MH5, MH7 and MH9 show- 
ing the most pronounced dip. During this phase, density 
perturbations seeded by the slightly supersonic turbulence 
grow more rapidly, and in these haloes the RMS density con- 
trast increases to well above unity. In MH5 and MH7, the 
density contrast even increases by an order of magnitude. 
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Figure 10. From left to right: number density of hydrogen nuclei, temperature, H2 fraction, and free-fall time over sound-crossing time 
in cubes of side length 200 au, centred on the densest cell in the halo and weighted with the mass and the square of the density of the 
cells traversed along the line of sight. From top to bottom: simulations employing 8, 16, 32, 64, and 128 cells per Jeans length according 
to [Turk et aT] ( |2010t in MHl. Except for an increasing amount of small-scale structure, there is no apparent trend in any of the shown 
quantities as the resolution is increased. 
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Figure 11. From top left to bottom right: number density of hydrogen nuclei, enclosed gas mass, temperature, H2 fraction, radial 
velocity over sound speed, rotation velocity over sound speed, rotation velocity over Keplerian velocity, turbulent Mach number, cooling 
time over free-fall time, effective equation of state, free-fall time over sound-crossing time, and RMS density contrast versus radius for 
8, 16, 32, 64, and 128 cells per Jeans length according to [Turk et al.| ( |2010| ) in MHl. The solid black line in the top left panel denotes 
an r~^-^ density profile. The grey dashed lines denote a cooling time equal to the free-fall time, an effective equation of state of 1.1, a 
free-fall time equal to the sound-crossing time, and an RMS density contrast of unity. Qualitatively, the profiles are similar among the 
various simulations, but also display a substantial amount of scatter. We find no clear trend in any of the profiles as the resolution is 
increased. 
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The formation of secondary clumps in these haloes is evi- 
dent from Figure [S] which shows the free-fall time over the 
sound-crossing time in the central 200 au of the box. Parcels 
of gas with a tendency to become gravitationally unstable 
are apparent. A comparison with Figure |6] shows that these 
regions are highly correlated with the isolated, dense clumps 
visible in this figure. 

Based on Figure [9] we investigate whether the clumps 
in MHl, MH5, MH7 and MH9 have accumulated enough 
mass to become gravitationally unstable, and whether they 
will collapse before they are accreted by the primary clump. 
To this end, the top panel shows the enclosed mass over 
the locally estimated Bonnor-Ebert mass (E bert|1955|[Boir 



1956), centred on the minimum in the free-fall time 



over the sound-crossing time as identified in Figure |8] The 
Bonnor-Ebert mass is computed as the mass- weighted aver- 
age among the cells within a given radius according to: 



1/2 /2^x3/2 



-3/2 2 



(12) 



The bottom panel shows the free-fall time over the accre- 
tion time. The free-fall time is computed by using the mass 
enclosed within a given radius, and the accretion time by 
tacc = d/vrad, whcrc d dcuotcs the distance of the clump 
from the center of the cloud (vertical lines in Figure |9|, and 
Vrad the net radial velocity of the mass enclosed within a 
given radius with respect to the center of the cloud. 

From this figure, it is evident that only the clouds 
in MH5 and MH7 are likely to fragment during the ini- 
tial free-fall phase. The enclosed mass exceeds the Bonnor- 
Ebert mass at small enough radii that the clumps are Jeans- 
unstable in their own right, which is not the case in MHl 
and MH9, where there is no clear distinction between the 
primary and secondary clumps. Furthermore, since the free- 
fall time falls below the accretion time at radii significantly 
smaller than the separation of the clumps from the center 
of the cloud, the clouds in MH5 and MH7 will likely col- 
lapse before they are accreted. However, the separation of 
the clumps is only of order a few tens of au - in contrast to 
the 800 au that were found in Turk et ah] ([2009 ) . Based on 
the results of Greif et al. (2012|, it thus appears likely that 



the secondary clumps will merge well before radiation feed- 
back can evacuate enough mass to alleviate the strong grav- 
itational torques that are present in the disc. On the other 
hand, the secondary clump found in [Turk et al.| ( |2012| has 
a much higher chance of surviving and forming a separate 
Population HI star. 



3.3 Resolution study 

Previous studies of Population HI star formation have ar- 
gued that the nature of primordial gas clouds depends sen- 
sitively on the employed resolution (e.g. Turk et aL]|2012 ). 
We investigate whether the resolution requirements found in 
these simulations also apply to the moving-mesh simulations 
presented here. To this end, we employ between 8 and 128 
cells per Jeans length according to Turk et al.^ (2010 ) , using 
the halo that was described in detail in Section 3.1 (MHl). 

In Figure [lO] we show the number density of hydro- 
gen nuclei, temperature, H2 fraction, and free-fall time over 
sound-crossing time in the central 200 au of the box. The 
morphology of the gas is highly irregular in all simulations, 




4 6 8 10 12 14 

log riH [cm"^] 

Figure 12. Temperature versus number density of hydrogen nu- 
clei for 8, 16, 32, 64, and 128 cells per Jeans length according 
to [Turk et aT] ( |2010| ) in MHl. The logarithm of the mass per bin 
over the total mass in the box is color-coded from blue (lowest) to 
red (highest), and the solid lines show the mass- weighted average 
value versus density. As the resolution is increased, the dispersion 
of the gas increases in particular in phases in which the gas cools 
efficiently, which is most likely due to the better sampling of the 
extreme ends of the distribution. 



and the amount of small-scale structure increases as the res- 
olution is increased. However, we do not find any other con- 
vincing trend in the density, temperature, H2 fraction, and 
free-fall time over sound-crossing time. Regions with a ten- 
dency to become gravitationally unstable are visible in all 
simulations, but their size or abundance does not decrease 
or increase as the resolution is increased. 

In Figure |11| we show radial profiles of the physical 
quantities discussed in Section 3.2. In the J 128 case, an ad- 
ditional resimulation was necessary to increase the dynamic 
range of the simulation, such that only the central 10^ au 
are shown. Qualitatively, the profiles agree well with each 
other, but there is also a substantial amount of scatter. As 
in Figure |10| we find no clear trend in any of the shown 
physical quantities as the resolution is increased. In Fig- 
ure [12] we show the distribution of the gas in density and 
temperature. For increasing resolution, the temperature dis- 
persion increases when the gas cools efficiently. This is most 
likely due to the better sampling of the extreme ends of the 
distribution, and may raise the probability for finding gravi- 
tationally unstable clumps of gas. However, we find no clear 
evidence for such a trend from Figures [lO] and [TT] Instead, 
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the random nature of self-gravitating collapse and the non- 
linear rate equations that govern the chemical and thermal 
evolution of the cloud appear to dominate over all other 
trends that may be present. In particular, the differences 
between simulations with different resolution are similar to 
the differences between simulations of different haloes. 



4 SUMMARY AND CONCLUSIONS 

We have investigated the operation of the chemo-thermal 
instability in primordial, star-forming clouds with a suite 
of nine high-resolution simulations. At densities nn ^ 
10^ cm"*^, three-body reactions rapidly increase the H2 frac- 
tion and trigger potential runaway cooling, which is counter- 
acted by the chemical heating of the gas due to the release 
of the binding energy of the H2 molecule, and the increasing 
optical depth to H2 line emission. The competition between 
these processes results in a drop of the cooling time below 
the free-fall time on a scale of ^ 1000 au, such that the ef- 
fective equation of state of the gas decreases from jef£ — 1.1 
to :^ 0.8. As a result, the free-fall time drops to well below 
the sound-crossing time on scales < 100 au, and perturba- 
tions seeded by the transonic turbulence within the cloud 
may grow large enough that they trigger the gravitational 
collapse of secondary clumps. In two of the nine haloes inves- 
tigated, enough mass accumulates in the secondary clumps 
that they will most likely collapse before they are accreted 
by the primary clump. 

In addition to the main simulation suite, we have per- 
formed a resolution study to investigate how our results de- 
pend on the number of cells employed per Jeans length. 
Apart from a broader temperature distribution of the gas, 
we find no convincing trend in the appearance of the cloud 
or other physical quantities as the resolution is increased. In- 
stead, the turbulent nature of self- gravitating collapse and 
the non-linear rate equations that govern the chemical and 
thermal evolution of the cloud appear to dominate over all 
other trends. The differences between simulations with dif- 
ferent resolution are comparable to the differences between 
simulations of different haloes. 

Qualitatively, our results agree well with previous three- 
dimensional simulations, but we also find a few differences. 



In the simulations of Yoshida et al. (2006), the primordial 



gas cloud remained marginally stable at high densities. Most 
haloes in our study follow this trend, but in some cases the 
transonic turbulence present in the gas seeds large enough 
perturbations that fragmentation may occur. Compared to 
Turk et al. (2009]), we find very good agreement concerning 
the relative fraction of haloes in which fragmentation oc- 
curs: in the present study, the gas fragments in two out of 
nine haloes, while in Turk et al. (2009) one out of five haloes 



fragments. However, this result may be somewhat mislead- 
ing, since fragmentation in our simulations occurs on a scale 
of a few tens of au, while in Turk et al. (2009) the secondary 
clump forms on a scale of ^ 1000 au. 

Apart from systematic differences concerning the grav- 
itational and hydrodynamic solvers, an additional source of 



discrepancy with respect to Turk et al. ( 2009 ) is the chem- 
ical model employed. The most important differences are 
the choice of the H2 cooling function, the three-body forma- 
tion rate, and a direction-dependent escape fraction for H2 



line emission, which raises the cooling efficiency of the gas 
( Turk et al.||2011 ). As a result, the temperature of the gas 
at 1 au from the center in our study is around 1000 K, 



while in |Turk et al.| ( |2009t it exceeds 2000 K. Collisional dis- 
sociation becomes effective before collision-induced emission 
cools the gas, and as a result the H2 fraction in [Turk et al.| 
(|2009 ) has begun to decrease at radii where we find that the 
gas is still fully molecular. Finally, Turk et al. ( ,2012J found 
that the gas becomes more rotationally supported when less 
than 32 cells per Jeans length are used. We do not find such 
a trend even when the resolution is decreased to only eight 
cells per Jeans length, which is close to the bare minimum 



required to resolve self- gravitating collapse (Truelove et al. 
1998||Greif et al.||201l| . At high resolution, both methods 
yield very similar cloud morphologies. 

One of the caveats of this study is that we have not 
taken supersonic streaming motions into account, which af- 
fects the mass scale and virialisation process of the haloes in 



which Population III star formation occurs ( 


Tseliakhovich 


& Hirata"2010| Greif et al.| 2011| |Maio et a 


I. 2011; Stacy 


et al.|2011i 


Naoz et al.| 2012| |2013| |0'Leary & McQuinn 


2012|). We also neglect the potential heating of the gas by 



self- annihilating DM ( 


Freese et al. 2008 locco et al. 2008 


Ripamonti et al. 2010 


Smith et al. 2012). Instead of per- 



forming full radiation transfer for H2 line emission, we resort 
to approximate, on-the-spot methods that likely reduce the 
cooling efficiency of the gas. Furthermore, some of the em- 



ployed chemical rate equations are highly uncertain ( Glover 
|fc Abel|[20Q8| [Turk et al.||2011|. We also neglect the influ- 
ence of magnetic fields, which might become dynamically 



important during the initial collapse ( 


Machida et al.||2008 


Xu et al.||2008|| 


Schleicher et al.||2009| |2010| Sur et al.||2010 


Federrath et al.|2011||Schober et al.|2012||Turk et al.|2012|). 



In concluding, we note that while only a fraction of 
the star-forming clouds in minihaloes are expected to frag- 
ment during their initial collapse, this by no means precludes 
fragmentation at a later stage, when the gas becomes rota- 
tionally supported in a Keplerian disc, and the infall rate 
is significantly reduced compared to simple one-dimensional 
estimates ( [Greif et al.||2012] ). The chemo-thermal instabil- 
ity and the resulting gravitational instability may therefore 
have substantially more time to develop than is apparent 
from the limited period of time simulated here. 
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